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Abstract 

Spatiotemporal properties of seismicity are investigated for a worldwide (WW) catalog and for South- 
ern California in the stationary case (SC), showing a nearly universal scaling behavior Distributions of 
distances between consecutive earthquakes (jumps) are magnitude independent and show two power-law 
regimes, separated by jump values about 200 km (WW) and 15 km (SC). Distributions of waiting times 
conditioned to the value of jumps show that both variables are correlated in general, but turn out to be in- 
dependent when only short or long jumps are considered. Finally, diffusion profiles reflect the shape of the 
jump distribution. 
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Earthquake statistics provides important insights about seismic processes 1 L Simple laws of 
seismicity such as the Gutenberg-Richter frequency-magnitude relation j^y, ^1 and the Omori 
law for the temporal decay of the seismic activity [4] allow for instance the prediction of the 
size of a characteristic earthquake, the evaluation of the spatial distribution of stress in a fault, 
or the estimation of seismic risk after a strong event ||5|, |g, LZD. It is natural to expect that more 
sophisticated, multidimensional statistical studies can be very valuable for hazard assessment as 
well as for understanding fundamental properties of earthquakes. In this way, the structure of 
earthquake occurrence in magnitude, space, and time should reflect all the complex interactions in 
the litho sphere. 

The first systematic analysis of seismicity taking into account its multidimensional nature was 
performed by Bak et al, who studied the dependence of waiting-time distributions on the magni- 
tude range and the size of the spatial region selected for analysis. Their study revealed, by means 
of a scaling law, the complexiW of a self-similar dynamical process with a hierarchy of structures 
over a wide range of scales isllsLuiilill]- Waiting times (also referred to as return times, inter- 
event times, etc.) are just the time intervals between consecutive earthquakes in a region, and can 
be studied in two ways: i) by using single-region waiting-time distributions, in which an arbitrary 
region is characterized by its own distribution [12], or ii) the original approach [8], in which a large 
region is divided into smaller, equally sized areas and waiting times are measured for the smaller 
areas but are included together into a unique mixed distribution for the whole region. (This pro- 
cedure can be interpreted as the measurement of first return times of seismic activity to different 
points in a large, spatially heterogeneous system, see Ref. [9].) The outcomes of both approaches 
are clearly different, but in any case scaling turns out to be a fundamental tool of analysis, reducing 
the multidimensional dependence of the waiting-time distributions to simple univariate functions. 
Further, it is possible to understand this scaling approach in terms of a renormalization-group 
transformation lisl . 

Equally important for risk estimations and forecasting, though much less studied ll^ fisl . 
should be the statistics of the distances between consecutive earthquakes, which we can identify 
with jumps (or flights) of earthquake occurrence. Very recently, Davidsen and Paczuski have 
provided a coherent picture using Bak et al.'s mixed-distributions procedure 1 16]; in contrast, our 
paper undertakes the study of the earthquake-distance problem considering the simpler approach 
of single-region distributions, and for the case of stationary seismicity. The results will lead us 
to examine the distribution of waiting times conditioned to different values of the jumps; finally. 
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we measure diffusion profiles for seismic occurrence, turning out to be that the profiles are largely 
determined by the distribution of jumps. 

Consider that an arbitrary spatial region and a range of magnitudes have been selected for 
analysis. The unit vector locating the epicenter of the z— th earthquake on the portion of the Earth 
surface selected is given by 

n = (cos 9/ cos 0/, sin (p,cos 0,-, sin Oi) , 

where the angles (pi and Oi denote longitude and latitude, respectively. The spatial distance, or 
jump, between the i—th event and the immediately previous in time event, z — 1, can be obtained 
from the angle a, defined by the two vectors, 

Ui = arccos(r,_i -r,). 

In this way one can measure distances as angles, in degrees; the distance in km is obtained by 
multiplying oc, in radians by the Earth radius (about 6370 km, assuming a spherical Earth). 

Given a set of values of the jumps, their probability density D{a) is defined as the probability 
per unit distance that the distance is in a small interval containing a. In order to avoid boundary 
problems it is convenient to start the analysis of seismicity on a global scale, which has also the 
advantage of stationarity (or, more properly, homogeneity in time, at least for the last 30 years). 
Stationarity means that any short time period is described by (roughly) the same seismic rate, 
R (defined as the number of earthquakes per unit time); in such a case, a linear increase of the 
cumulative number of earthquakes versus time must be observed. Notice that stationarity does 
not mean that aftershock sequences are not present in the data, rather, many sequences can be 
interwinned but without a predominant one in the spatial scale of observation. 

n 

The results for the NEIC-PDE worldwide catalog II17I1 . covering the period 1973-2002, are 
shown in Fig. [ija) (bottom set of curves), using events with magnitude M above different thresh- 
olds Mc (i.e., M > Mc). An important conclusion drawn from the figure is the independence of the 
jump distributions on the magnitude threshold, as in Ref. [ 16], implying that the spatial occurrence 
of large earthquakes is no different than the occurrence of small ones, in contrast to claims by some 
authors lisi . There is, nevertheless, an exception for small distances, less than about 0.1° (~ 10 
km) for Mc from 5.5 to 6, comparable with the size of rupture lioi . The distributions are charac- 
terized by a decreasing power-law regime from about 0.1° to 2°, with an exponent about 1.6, and 
a possible increasing second power law for a > 2°, with exponent 0.3, decaying abruptly close to 
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the maximum distance. Note the difference between these exponents and the value 0.6, obtained 
in Ref. 

mixing many areas in Southern California and including non- stationary periods. 
However, the measurement of D{a) contains an inconsistency when it is used for distances 
over a sphere; for instance, after any event, there are many ways to obtain a jump of a = 90° 
(360 "ways" in intervals of 1°) but there is only one way to get a = 180° (which is the maximum 
possible distance, corresponding to the antipodes of the initiating event). Therefore, instead of 
working with the probability density defined per unit distance, i.e., per unit angle, it is preferable 
to use (and easier to interpret) the probability density defined per unit solid angle (taking into 
account all the possible orientations for a given distance a). A straightforward way to estimate 
this quantity is by means of the transformation 

D{a)^D{a)=D{a)/{2nsma), 

due to the fact that iKs'mada is the element of solid angle defined by the points at distance a. 

The results of the transformation are clearly seen in Fig. [Ha) (top curves). In addition to the 
power law for the range 0.1° — 2°, whose exponent turns out to be ~ 2.6, the second power law, 
now decreasing, becomes more apparent from about 2° to the maximum distance, 180°, with an 
exponent ~ 0.7. This second power law seems to imply a dependence of events at large spatial 
scales. Indeed, for a given direction, it is more likely a jump with a = 20° than one of 160° 
(independence would imply that the probabilities were the same); however, this effect is due to 
the fact that earthquake occurrence is not uniform over the Earth but fractal, and the power law 
reflects the fractal structure of the epicenters. If for long distances there is no causal relation 
between events, the distribution of earthquake jumps is equivalent to the distribution of distances 
between any pair of earthquakes (not necessarily consecutive); and as the density corresponding 
to this distance (per unit distance) is the derivative of the well-known correlation integral, widely 
used to obtain fractal dimensions, independence at long distances implies 

D(a) ~ a''f'^ and D{a) ~ 1 /a^^''f for long a, 

where df is the correlation dimension. We have measured the probability density of the distance 
between any two earthquakes (per unit solid angle), obtaining a behavior proportional to 1/ a^-^^, 
for events with M > 6, which implies a correlation dimension <i/ — 1 . 15 in reasonable agreement 
with D{a) ~ 1/ a^-^, and confirming the independence of events for a > 2°. 

Turning back to the short-jump regime, the excess of probability given by the corresponding 
power-law with respect the long -jump power law (associated to the uncorrelated regime) is a clear 
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sign of spatial clustering in earthquake occurrence, extending up to distances of about 2° ~ 200 
km. It is likely that this clustering extends beyond this limit, but it is not detectable with this 
procedure as it is hidden below the uncorrelated, long-distance regime. In any case, the power law 
behavior implies the no existence of a finite correlation length, at least up to 200 km, at variance 



with the findings of Ref. ||2C , 



2111 . In consequence we can consider the earthquake jumps as Levy 
flights, which are one of the main causes of anomalous diffusion. 

In contrast to worldwide seismicity, regional seismicity turns out to be nonstationary, in gen- 
eral, and in the same way as for waiting-time distributions 13], the distributions of earthquake 
jumps depend on the time window selected for analysis. This problem will be avoided here by 
considering specific time windows characterized by stationary seismicity; an important realization 
then is that stationary seismicity is characterized by stationary distributions. Additionally, sta- 
tionary seismicity has the advantage of a lower magnitude threshold of completeness, due to the 
absence of large events. 

We consider Southern-California seismicity, from the waveform cross-correlation catalog by 
Shearer et al. I22I1 . The analysis has been performed on 1 1 stationary periods, comprised between 
the years 1986 and 2002, yielding a total time span of 9.25 years and containing 6072 events for 
M > 2.5. It is very striking, as Fig. [Tib) shows, that the distributions of jumps resemble very much 
those of the worldwide case, with two power-law regimes, but in a different scale (the separation 
is at about 0.1°). In principle, for a smaller area such as California one would expect that the 
distribution of jumps is just a truncation of the global one, but Fig. [Tib) clearly refutes this fact, 
providing a clear illustration of earthquake-occurrence self-similarity in space. In fact, the figure 
shows these distributions rescaled by a factor L, with L the maximum distance for the region, 
which is L ~ 6.5° for Southern California, and L = 180° for the worldwide case (also included in 
the plot). The reasonable data collapse is a signature of the existence of a scaling law for the jump 
distribution, 

D{a)-g{a/L)/L 

(for D{a) we need to include an extra factor C of normalization). The scaling law can only be 
given the status of approximated, as the exponent for short jumps for Southern California seems 
to be different than in the worldwide case, 2.15 versus 2.6 for the distribution D. Nevertheless, 
this variation could be due to artifacts in the short-distance properties of the catalogs. On the 
other hand, the long-jump regime is well described by the worldwide exponent. We recall that 
the exponents are far from the value of Davidsen and Paczuski for the nonstationary California 



case 1 161. The main difference with the worldwide case is the abrupt decay of D{a) close to the 
maximum a. Clearly, this is a boundary effect. We have tried several corrections for this effect (in 
particular the factor 2;rsinci5 is not right when the ring defined by the set of points at a distance a 
is not totally included in the region under study) but we have not found a substantial improvement 
and then we have preferred to keep the things simple. 

The distributions of jumps we have obtained, together with previous work on the distribution 
of waiting times 13], provide a simple picture of earthquake occurrence as a continuous-time 
random walk. This means that we can understand seismic activity as a random walk over the Earth 
surface, where one event comes after another at a distance that follows D{a) and after a waiting 
time T given by the waiting-time probability density, D{t). However, an important ingredient is 
missing in this picture, as we need to take into account the correlations between distances and 
waiting times. 

With this goal in mind, we introduce the conditional waiting-time probability density, D{t\X), 
where | X means that only the cases where X is verified are taken into account; in our case X will 
be a set of values of the jumps. It turns out that if Z)(t | X) does not depend on X, then T and X are 
independent, whereas when D{t\X) changes with X, then T and X are correlated (nonlinearly, in 
general). 

From the behavior of D{a), a natural threshold for the jumps in the worldwide case is a ~ 2°. 
Figure |2ta) compares D{x\a short ) with D{x\a long ) for worldwide seismicity, where short 
and long refer to sets of distances below and above 2°, respectively. The differences are clear: for 
short jumps the waiting-time distribution is a decreasing power law for several decades, with an 
exponent close to 1, and ends in an exponential decay. We can identify these distributions with 
highly correlated events, i.e., aftershock sequences. On the other hand, for long a the waiting- 
time distribution seems to be exponential for its full range, compatible with a Poisson process and 
therefore with independent occurrence. 

But more surprising than the differences between short and long jumps are perhaps the similar- 
ities for short jumps. In fact, there seems to be a radical change of behavior separated by a ~ 2°, 
in the sense that if we are above, or below, this value, the distributions do not change, in other 
words, D{x\a < 0.25°) = D(t|0.25 <a< 1°), etc. and D(t|2° < a < 30°) = D{x\a > 100°), 
etc., see Fig. |2ta) (for the range between 1° and 2° the behavior is not clear as the statistics is 
low). This means that for short jumps the waiting-time distribution is independent on the value of 
the jump, and the same happens for long jumps, but when the whole range of jumps is considered 
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this is no longer true and both variables become dependent, in contrast to Ref. 11^ . For each set 
of curves we can fit a gamma distribution, D{T\a) oc e^^/'^/r^^^, turning out to be 7 ~ 0.17 for 
short a and 7^ 0.9 for long a. The latter value is very close to 1, the characteristic value of a 
Poisson process; in fact, the difference between 1 and 0.9 is not significant within the uncertainty 
of our data, but if the hypothesis of a Poisson process for long distances could be rejected at any 
reasonable significance level (for which much more data would be necessary) and a value 7 < 1 
could be significantly established, this would constitute a support for the existence of long-range 
earthquake triggering. In any case, our findings are in disagreement with the hypothesis put for- 
ward in Ref. Jioll . which considers seismicity as a process uncorrelated in time but correlated in 
space. 

This result confirms that the universal scaling law found for stationary seismicity in Ref. I12I1 
is in fact a mixture of aftershocks and independent events, turning out to be very striking that this 
mixture leads to a universal behavior. The reason could be a universal proportion of aftershocks 
versus mainshocks in stationary seismicity lisl . Further, the conditional time distributions verify 
a scaling law, 

D(T|a±) =i?(^,M,)/±(i?(=^,M,)T), 

where the indices -|- and — refer to long and short jumps, respectively, and R^^{M,Mc) is the mean 
waiting time of the unconditional distribution for the region ^ and M > Mc (obviously, a scaling 
with the mean of the conditional distribution also holds). FigureEtb) shows this behavior. 

The next step is to measure earthquake diffusion directly from data. The fundamental quantity 
is p(aj), which we call diffusion profile and gives the probability density that two earthquakes 
(not consecutive) are at a distance a when they are separated by a time t, i.e., 

p{a,t)da = Prob[a < distance at a time t < a + da]. 

In practice, the single value of t is replaced by a range of values and, as above, we introduce 
p{a,t) = p{a,t)/{2ns'ma). In the case of normal diffusion p{a,t) is given by a Gaussian dis- 
tribution (more precisely, a semi-Gaussian, as a > 0), with a second moment scaling as (oc^) ~ t, 
whereas p(a, t) is given by a Rayleigh distribution (or a Maxwell distribution in three dimensions). 
In contrast, our measurements yield to results far from normality, see Fig. |3l For short times, the 
profile resembles the distribution of jumps with 2 regimes; as time evolves events migrate farther 
from the origin, towards the long-range part of the curve, which becomes dominant for all a for 
long times, with an exponent of value 0.8 for p. 
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A scaling law also holds for the diffusion profile, i.e., 



p{a,t) = h{a/L,t)/L. 



Finally, the process turns out to be clearly subdiffusive, as (a^) grows more s 



owly than lin- 



early, although we cannot provide a definite value of a exponent as in Refs. 11241 

This paper would not have been possible without the effort of the researchers who have made 
earthquake catalogs publicly available. The Ramon y Cajal program and the projects BFM2003- 
06033 (MCyT) and 2001-SGR-00186 (DGRGC) are acknowledged. 
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FIG. 1: (color online) (a) Probability densities of earthquake jumps, defined per unit distance (angle), 
D{a), and per unit solid angle, D{a) = D(a)/(360°sina) (shifted up a factor 1000 for clarity sake), for 
worldwide seismicity with lower magnitudes from M^- = 5 to = 6.5. Lines come from power-law fits to 
D(a) for Mc = 5. (b) Rescaled probability densities D{a) for worldwide seismicity (WW) and for several 
stationary periods in Southern California (SC), with different Mc values. For each case, L = 180° and 
L = 6.5°, and the normalization factors C = 0.0235 and C = 1.14 degrees^' are determined for Mc = 5 and 
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FIG. 2: (color online) (a) Waiting-time probability densities conditioned to different sets of values of the 
jumps, in particular a short (below 1°) and a long (above 2°), using worldwide seismicity withM > 5.5. 
In each case, gamma fits to all the curves are shown, with parameter 7 = 0. 17 (a short) and 7 = 0.90 (a 
long), (b) Rescaled waiting-time probability densities conditioned to short jumps, a < 1° for worldwide 
seismicity (WW) and a < 0. 1° for Southern California (SC). Different Mc values are used. The rescaling 
factor R refers to the unconditional distribution, and only depends on Mc and on the spatial region The 
solid line is the gamma fit in (a) rescaled in the same way. 
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FIG. 3: (color online) Diffusion profile p{(X,t) versus a (rescaled by L) for worldwide seismicity with 
M > 5.5 (bottom curves) and for Southern California with M > 2.5 (top). Time ranges from tens of sec- 
onds to several years. The displayed power laws are indicative, and have the same exponents as the jump 
distributions. 
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